{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from linearmodels.panel import PanelOLS\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import LabelEncoder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Analysis 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1. Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"C:\\\\df_analysis_1.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assuming df is your DataFrame and it's already loaded\n",
    "df = df.set_index(['si_gun_gu', 'year'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to run fixed effects regression\n",
    "def run_fe_regression(data, dependent_var, independent_vars):\n",
    "    y = data[dependent_var]\n",
    "    X = data[independent_vars]\n",
    "    X = sm.add_constant(X)\n",
    "    model = PanelOLS(y, X, entity_effects=True)\n",
    "    results = model.fit(cov_type='clustered', cluster_entity=True)\n",
    "    return results\n",
    "\n",
    "# Define independent variables\n",
    "independent_vars = ['Urbanization', 'Income', 'Homeownership', 'Air_pollution']\n",
    "\n",
    "# Models for Distrust\n",
    "results_distrust_evi = run_fe_regression(df, 'Distrust_1', ['EVI'] + independent_vars)\n",
    "results_distrust_ndvi = run_fe_regression(df, 'Distrust_1', ['NDVI'] + independent_vars)\n",
    "\n",
    "# Models for Exploitation\n",
    "results_exploitation_evi = run_fe_regression(df, 'Exploitation_1', ['EVI'] + independent_vars)\n",
    "results_exploitation_ndvi = run_fe_regression(df, 'Exploitation_1', ['NDVI'] + independent_vars)\n",
    "\n",
    "# Models for Indifference\n",
    "results_indifference_evi = run_fe_regression(df, 'Indifference_1', ['EVI'] + independent_vars)\n",
    "results_indifference_ndvi = run_fe_regression(df, 'Indifference_1', ['NDVI'] + independent_vars)\n",
    "\n",
    "\n",
    "# Collect results into a list\n",
    "results_list = [\n",
    "    results_distrust_evi, results_distrust_ndvi, \n",
    "    results_exploitation_evi, results_exploitation_ndvi,\n",
    "    results_indifference_evi, results_indifference_ndvi\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract coefficients and CIs\n",
    "def extract_coefs(results, var):\n",
    "    coefs = results.params[var]\n",
    "    ci = results.conf_int().loc[var]\n",
    "    ci_lower = ci[0]\n",
    "    ci_upper = ci[1]\n",
    "    return coefs, ci_lower, ci_upper\n",
    "\n",
    "# Define variables and models\n",
    "variables = ['EVI', 'NDVI']\n",
    "models = ['Distrust', 'Exploitation', 'Indifference']\n",
    "\n",
    "# Extract coefficients and CIs\n",
    "data_for_plot = []\n",
    "for i, var in enumerate(variables):\n",
    "    for j, model in enumerate(models):\n",
    "        result = results_list[j*2 + i]  # Corrected indexing\n",
    "        coefs, ci_lower, ci_upper = extract_coefs(result, var)\n",
    "        data_for_plot.append([var, f'Model {model}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "# Create DataFrame\n",
    "df_plot = pd.DataFrame(data_for_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(12, 6))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Variables for adjusting positions\n",
    "offset = 0.2\n",
    "positions = {'EVI': -offset, 'NDVI': 0,}\n",
    "colors = {'EVI': 'blue', 'NDVI': 'green', }\n",
    "\n",
    "# Plotting each variable\n",
    "for var in variables:\n",
    "    subset = df_plot[df_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    adjusted_positions = [x + positions[var] for x in x_positions]\n",
    "    \n",
    "    plt.errorbar(adjusted_positions, subset['Coef'], yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                 fmt='o', label=var, capsize=5, markersize=18, color=colors[var])\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker\n",
    "    for j in range(len(subset)):\n",
    "        plt.text(adjusted_positions[j], subset['Coef'].iloc[j] + 0.5, f'{subset[\"Coef\"].iloc[j]:.3f}', \n",
    "                 ha='center', va='bottom', fontsize=24, color=colors[var])\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "plt.axhline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['Model 1\\n(DV: Distrust)', 'Model 2\\n(DV: Exploitation)', 'Model 3\\n(DV: Indifference)']\n",
    "plt.xticks(range(len(custom_labels)), custom_labels, size=22)\n",
    "plt.yticks(size=16)\n",
    "\n",
    "# plt.xlabel('Model')\n",
    "# plt.ylabel('Coefficient')\n",
    "# plt.title('Coefficient Plot')\n",
    "\n",
    "# Increasing legend text size\n",
    "plt.legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large')\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2. Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"C:\\\\df_analysis_1.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assuming df is your DataFrame and it's already loaded\n",
    "df = df.set_index(['si_gun_gu', 'year'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to run fixed effects regression\n",
    "def run_fe_regression(data, dependent_var, independent_vars):\n",
    "    y = data[dependent_var]\n",
    "    X = data[independent_vars]\n",
    "    X = sm.add_constant(X)\n",
    "    model = PanelOLS(y, X, entity_effects=True)\n",
    "    results = model.fit(cov_type='clustered', cluster_entity=True)\n",
    "    return results\n",
    "\n",
    "# Define independent variables\n",
    "independent_vars = ['Urbanization', 'Income', 'Homeownership', 'Air_pollution']\n",
    "\n",
    "# Models for Distrust\n",
    "results_distrust_evi = run_fe_regression(df, 'Distrust_1', ['EVI'] + independent_vars)\n",
    "results_distrust_ndvi = run_fe_regression(df, 'Distrust_1', ['NDVI'] + independent_vars)\n",
    "\n",
    "# Models for Exploitation\n",
    "results_exploitation_evi = run_fe_regression(df, 'Exploitation_1', ['EVI'] + independent_vars)\n",
    "results_exploitation_ndvi = run_fe_regression(df, 'Exploitation_1', ['NDVI'] + independent_vars)\n",
    "\n",
    "# Models for Indifference\n",
    "results_indifference_evi = run_fe_regression(df, 'Indifference_1', ['EVI'] + independent_vars)\n",
    "results_indifference_ndvi = run_fe_regression(df, 'Indifference_1', ['NDVI'] + independent_vars)\n",
    "\n",
    "# Generate predicted effects\n",
    "evi_range = np.linspace(0.1, 0.4, 100)\n",
    "ndvi_range = np.linspace(0.2, 0.7, 100)\n",
    "\n",
    "def predict_effects(results, var_name, var_range, independent_vars):\n",
    "    predictions = []\n",
    "    for value in var_range:\n",
    "        data = {var_name: value, 'const': 1}\n",
    "        for var in independent_vars:\n",
    "            data[var] = df[var].mean()  # Using mean of other independent vars\n",
    "        X_new = pd.DataFrame(data, index=[0] * len(df.index))\n",
    "        X_new.index = df.index  # Set the same index as the original data\n",
    "        X_new = X_new[['const'] + [var_name] + independent_vars]\n",
    "        pred = results.predict(exog=X_new)\n",
    "        predictions.append(pred.mean().item())  # Use mean to aggregate the predictions\n",
    "    return predictions\n",
    "\n",
    "predictions_distrust_evi = predict_effects(results_distrust_evi, 'EVI', evi_range, independent_vars)\n",
    "predictions_distrust_ndvi = predict_effects(results_distrust_ndvi, 'NDVI', ndvi_range, independent_vars)\n",
    "\n",
    "predictions_exploitation_evi = predict_effects(results_exploitation_evi, 'EVI', evi_range, independent_vars)\n",
    "predictions_exploitation_ndvi = predict_effects(results_exploitation_ndvi, 'NDVI', ndvi_range, independent_vars)\n",
    "\n",
    "predictions_indifference_evi = predict_effects(results_indifference_evi, 'EVI', evi_range, independent_vars)\n",
    "predictions_indifference_ndvi = predict_effects(results_indifference_ndvi, 'NDVI', ndvi_range, independent_vars)\n",
    "\n",
    "# Plotting the results\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "\n",
    "# Plot for EVI\n",
    "axes[0].plot(evi_range, predictions_distrust_evi, label='Distrust', color='blue', linestyle='-', linewidth=2)\n",
    "axes[0].plot(evi_range, predictions_exploitation_evi, label='Exploitation', color='green', linestyle='--', linewidth=2)\n",
    "axes[0].plot(evi_range, predictions_indifference_evi, label='Indifference', color='red', linestyle=':', linewidth=2)\n",
    "# Annotate only the first and last points\n",
    "axes[0].text(evi_range[0], predictions_distrust_evi[0], f'{predictions_distrust_evi[0]:.2f}', fontsize=15, color='blue', ha='right', va='bottom')\n",
    "axes[0].text(evi_range[-1], predictions_distrust_evi[-1]+0.07, f'{predictions_distrust_evi[-1]:.2f}', fontsize=15, color='blue', ha='left', va='top')\n",
    "axes[0].text(evi_range[0], predictions_exploitation_evi[0]-0.07, f'{predictions_exploitation_evi[0]:.2f}', fontsize=15, color='green', ha='right', va='bottom')\n",
    "axes[0].text(evi_range[-1], predictions_exploitation_evi[-1]-0.04, f'{predictions_exploitation_evi[-1]:.2f}', fontsize=15, color='green', ha='left', va='top')\n",
    "axes[0].text(evi_range[0], predictions_indifference_evi[0]+0.08, f'{predictions_indifference_evi[0]:.2f}', fontsize=15, color='red', ha='right', va='bottom')\n",
    "axes[0].text(evi_range[-1], predictions_indifference_evi[-1]-0.16, f'{predictions_indifference_evi[-1]:.2f}', fontsize=15, color='red', ha='left', va='top')\n",
    "axes[0].set_xlabel('EVI', fontsize=16)\n",
    "axes[0].set_ylabel('Predicted Value', fontsize=14)\n",
    "axes[0].legend(fontsize=15)\n",
    "axes[0].set_title('(a) Predicted Effects for EVI', fontsize=18, pad=22)\n",
    "\n",
    "# Plot for NDVI\n",
    "axes[1].plot(ndvi_range, predictions_distrust_ndvi, label='Distrust', color='blue', linestyle='-', linewidth=2)\n",
    "axes[1].plot(ndvi_range, predictions_exploitation_ndvi, label='Exploitation', color='green', linestyle='--', linewidth=2)\n",
    "axes[1].plot(ndvi_range, predictions_indifference_ndvi, label='Indifference', color='red', linestyle=':', linewidth=2)\n",
    "# Annotate only the first and last points\n",
    "axes[1].text(ndvi_range[0], predictions_distrust_ndvi[0]-0.2, f'{predictions_distrust_ndvi[0]:.2f}', fontsize=15, color='blue', ha='right', va='bottom')\n",
    "axes[1].text(ndvi_range[-1], predictions_distrust_ndvi[-1]+0.08, f'{predictions_distrust_ndvi[-1]:.2f}', fontsize=15, color='blue', ha='left', va='top')\n",
    "axes[1].text(ndvi_range[0], predictions_exploitation_ndvi[0]+0.07, f'{predictions_exploitation_ndvi[0]:.2f}', fontsize=15, color='green', ha='right', va='bottom')\n",
    "axes[1].text(ndvi_range[-1], predictions_exploitation_ndvi[-1], f'{predictions_exploitation_ndvi[-1]:.2f}', fontsize=15, color='green', ha='left', va='top')\n",
    "axes[1].text(ndvi_range[0], predictions_indifference_ndvi[0], f'{predictions_indifference_ndvi[0]:.2f}', fontsize=15, color='red', ha='right', va='bottom')\n",
    "axes[1].text(ndvi_range[-1], predictions_indifference_ndvi[-1], f'{predictions_indifference_ndvi[-1]:.2f}', fontsize=15, color='red', ha='left', va='top')\n",
    "axes[1].set_xlabel('NDVI', fontsize=16)\n",
    "axes[1].set_ylabel('Predicted Value', fontsize=14)\n",
    "axes[1].legend(fontsize=15)\n",
    "axes[1].set_title('(b) Predicted Effects for NDVI', fontsize=18, pad=22)\n",
    "\n",
    "\n",
    "# Increase the size of x-ticks and y-ticks\n",
    "for ax in axes:\n",
    "    ax.tick_params(axis='x', labelsize=14)\n",
    "    ax.tick_params(axis='y', labelsize=14)\n",
    "\n",
    "# Adjust layout for better spacing\n",
    "plt.subplots_adjust(wspace=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Analysis 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"C:\\Users\\Han\\Desktop\\document\\불평등 연구\\주제36_녹지불평등\\npj urban sustainability\\df_merge.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop rows with missing values for simplicity\n",
    "df = df.dropna(subset=['Loneliness', 'EVI', 'NDVI', 'sex', 'marital', 'age', 'education', 'income', 'homeownership', 'religion', 'year'])\n",
    "\n",
    "# Ensure proper encoding of categorical variables\n",
    "df['sido_1'] = LabelEncoder().fit_transform(df['sido'])\n",
    "df['si_gun_gu_1'] = LabelEncoder().fit_transform(df['si_gun_gu'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify the grouping variable\n",
    "assert df['si_gun_gu_1'].isnull().sum() == 0, \"Grouping variable has missing values.\"\n",
    "assert len(df['si_gun_gu_1'].unique()) > 1, \"Grouping variable should have more than one unique value.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# Function to run mixed effects model\n",
    "def run_mixed_model(data, formula, re_formula):\n",
    "    model = smf.mixedlm(formula, data, groups=data[\"si_gun_gu_1\"], re_formula=re_formula)\n",
    "    result = model.fit()\n",
    "    return result\n",
    "\n",
    "# Mixed models\n",
    "formulas = [\n",
    "    'Loneliness ~ EVI + sex + marital + age + education + income + homeownership + religion + C(year)',\n",
    "    'Loneliness ~ NDVI + sex + marital + age + education + income + homeownership + religion + C(year)'\n",
    "]\n",
    "re_formulas = ['~EVI', '~NDVI']\n",
    "\n",
    "# Run mixed models\n",
    "results_loneliness = [run_mixed_model(df, formulas[i], re_formulas[i]) for i in range(2)]\n",
    "\n",
    "# Display the results\n",
    "for i, result in enumerate(results_loneliness):\n",
    "    print(f\"Model {i+1}\")\n",
    "    print(result.summary())\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to run mixed effects model\n",
    "def run_mixed_model(data, formula, re_formula):\n",
    "    model = smf.mixedlm(formula, data, groups=data[\"si_gun_gu_1\"], re_formula=re_formula)\n",
    "    result = model.fit()\n",
    "    return result\n",
    "\n",
    "# Mixed models\n",
    "formulas = [\n",
    "    'Loneliness ~ EVI + education + income + homeownership + C(year)',\n",
    "    'Loneliness ~ NDVI + education + income + homeownership + C(year)'\n",
    "]\n",
    "re_formulas = ['~EVI', '~NDVI']\n",
    "\n",
    "# Run mixed models\n",
    "results_loneliness = [run_mixed_model(df, formulas[i], re_formulas[i]) for i in range(2)]\n",
    "\n",
    "# Display the results\n",
    "for i, result in enumerate(results_loneliness):\n",
    "    print(f\"Model {i+1}\")\n",
    "    print(result.summary())\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract coefficients and CIs\n",
    "def extract_coefs(results, var):\n",
    "    coefs = results.params[var]\n",
    "    ci = results.conf_int().loc[var]\n",
    "    ci_lower = ci[0]\n",
    "    ci_upper = ci[1]\n",
    "    return coefs, ci_lower, ci_upper\n",
    "\n",
    "# Extract for each model\n",
    "variables = ['EVI', 'NDVI']\n",
    "data_for_plot = []\n",
    "for i, var in enumerate(variables):\n",
    "    coefs, ci_lower, ci_upper = extract_coefs(results_loneliness[i], var)\n",
    "    data_for_plot.append([var, f'Model {i+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_plot = pd.DataFrame(data_for_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mixed models for Distrust, Exploitation, and Indifference\n",
    "additional_formulas = [\n",
    "    'Distrust ~ Loneliness + sex + marital + age + education + income + homeownership + religion + C(year)',\n",
    "    'Exploitation ~ Loneliness + sex + marital + age + education + income + homeownership + religion + C(year)',\n",
    "    'Indifference ~ Loneliness + sex + marital + age + education + income + homeownership + religion + C(year)'\n",
    "]\n",
    "\n",
    "results_additional = [run_mixed_model(df, formula, '~NDVI') for formula in additional_formulas]\n",
    "\n",
    "# Display the results\n",
    "for i, result in enumerate(results_additional):\n",
    "    print(f\"Model {i+1}\")\n",
    "    print(result.summary())\n",
    "    print(\"\\n\")\n",
    "\n",
    "# Extract and plot coefficients for these models as done before.\n",
    "data_for_additional_plot = []\n",
    "for i, var in enumerate(['Loneliness']*3):\n",
    "    coefs, ci_lower, ci_upper = extract_coefs(results_additional[i], var)\n",
    "    data_for_additional_plot.append([var, f'Model {i+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_additional_plot = pd.DataFrame(data_for_additional_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mixed models for Distrust, Exploitation, and Indifference\n",
    "additional_formulas = [\n",
    "    'Distrust ~ Loneliness + education + income + homeownership + C(year)',\n",
    "    'Exploitation ~ Loneliness + education + income + homeownership + C(year)',\n",
    "    'Indifference ~ Loneliness + education + income + homeownership + C(year)'\n",
    "]\n",
    "\n",
    "results_additional = [run_mixed_model(df, formula, '~NDVI') for formula in additional_formulas]\n",
    "\n",
    "# Display the results\n",
    "for i, result in enumerate(results_additional):\n",
    "    print(f\"Model {i+1}\")\n",
    "    print(result.summary())\n",
    "    print(\"\\n\")\n",
    "\n",
    "# Extract and plot coefficients for these models as done before.\n",
    "data_for_additional_plot = []\n",
    "for i, var in enumerate(['Loneliness']*3):\n",
    "    coefs, ci_lower, ci_upper = extract_coefs(results_additional[i], var)\n",
    "    data_for_additional_plot.append([var, f'Model {i+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_additional_plot = pd.DataFrame(data_for_additional_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assuming df_plot and df_additional_plot are already defined with the necessary data\n",
    "\n",
    "# Define the plot size and style\n",
    "fig, axes = plt.subplots(1, 2, figsize=(24, 10))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Variables for adjusting positions\n",
    "offset = 0.3\n",
    "variables = ['EVI', 'NDVI']  # Define the variables to be plotted\n",
    "positions = {'EVI': -offset, 'NDVI': +offset}\n",
    "colors = {'EVI': 'blue', 'NDVI': 'green'}\n",
    "\n",
    "# Adjust the offset for annotations to avoid overlap\n",
    "annotation_offsets = {'EVI': 0.1, 'NDVI': -0.1}\n",
    "\n",
    "# Plot 1: Loneliness with EVI, NDVI\n",
    "for var in variables:\n",
    "    subset = df_plot[df_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    adjusted_positions = [x + positions[var] for x in x_positions]\n",
    "    \n",
    "    axes[0].errorbar(adjusted_positions, subset['Coef'], \n",
    "                     yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                     fmt='o', label=var, capsize=5, markersize=22, color=colors[var])\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker, with offset to avoid overlap\n",
    "    for j in range(len(subset)):\n",
    "        axes[0].text(adjusted_positions[j], \n",
    "                     subset['Coef'].iloc[j] + 3.5 + annotation_offsets[var],  # Add offset to separate annotations\n",
    "                     f'{subset[\"Coef\"].iloc[j]:.3f}', ha='center', va='bottom', \n",
    "                     fontsize=26, color=colors[var])\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "axes[0].axhline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "axes[0].set_xlim(-1, 1)  # Adjust x-axis limits to center the label\n",
    "axes[0].set_xticks([0])  # Position tick in the center\n",
    "axes[0].set_xticklabels(['DV: Loneliness'], size=32, ha='center')\n",
    "axes[0].tick_params(axis='y', labelsize=22)\n",
    "\n",
    "# Increasing legend text size\n",
    "axes[0].legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large')\n",
    "axes[0].grid(True)\n",
    "\n",
    "# Title for left plot\n",
    "axes[0].set_title('(a) Vegetation Index and Loneliness', fontsize=41, pad=29)\n",
    "\n",
    "# Plot 2: Loneliness with NDVI for Distrust, Exploitation, and Indifference\n",
    "for var in ['Loneliness']:\n",
    "    subset = df_additional_plot[df_additional_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    \n",
    "    axes[1].errorbar(x_positions, subset['Coef'], \n",
    "                     yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                     fmt='o', label=var, capsize=5, markersize=22, color='red')\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker\n",
    "    for j in range(len(subset)):\n",
    "        axes[1].text(x_positions[j], subset['Coef'].iloc[j] + 0.004, \n",
    "                     f'{subset[\"Coef\"].iloc[j]:.3f}', ha='center', va='bottom', \n",
    "                     fontsize=26, color='red')\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "axes[1].axhline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['Model 1\\n(DV: Distrust)', 'Model 2\\n(DV: Exploitation)', 'Model 3\\n(DV: Indifference)']\n",
    "axes[1].set_xticks(range(len(custom_labels)))\n",
    "axes[1].set_xticklabels(custom_labels, size=28)\n",
    "axes[1].tick_params(axis='y', labelsize=22)\n",
    "\n",
    "# Increasing legend text size and moving it to the bottom right\n",
    "axes[1].legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large', loc='lower right', bbox_to_anchor=(1, 0))\n",
    "axes[1].grid(True)\n",
    "\n",
    "# Title for right plot\n",
    "axes[1].set_title('(b) Loneliness and Fragmentation', fontsize=41, pad=29)\n",
    "\n",
    "plt.tight_layout(pad=5.0)  # Increase space between plots\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
